Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update documentation #16

Merged
merged 18 commits into from
Sep 13, 2023
Merged

Update documentation #16

merged 18 commits into from
Sep 13, 2023

Conversation

FredericWantiez
Copy link
Member

No description provided.

@FredericWantiez FredericWantiez changed the title Fred/doc Update documentation Jul 6, 2023
examples/smc.jl Outdated Show resolved Hide resolved
examples/smc.jl Outdated Show resolved Hide resolved
src/SSMProblems.jl Outdated Show resolved Hide resolved
@FredericWantiez
Copy link
Member Author

573cab3 actually feels a bit weird, the model definition is now specialized to the concrete particle type:

function transition!!(rng::AbstractRNG, t::Int, particle::Particle{<:LinearSSM})
    if t == 1
        return Particle(LinearSSM(rand(rng, f0(t))))
    else
        return Particle(LinearSSM(rand(rng, f(t, particle.model.state))))
    end
end

where if we just dispatch on the model, like in 6394704, we can drop-in different particle estimators with different specialization of AbastractParticle{T}

@yebai
Copy link
Member

yebai commented Jul 28, 2023

function transition!!(rng::AbstractRNG, t::Int, particle::Particle{<:LinearSSM})
    if t == 1
        return Particle(LinearSSM(rand(rng, f0(t))))
    else
        return Particle(LinearSSM(rand(rng, f(t, particle.model.state))))
    end
end


function transition!!(rng::AbstractRNG, t::Int, m::LinearSSM, particle::Particle)
    if t == 1
        return Particle(rand(rng, f0(t)), m)
    else
        return Particle(rand(rng, f(t, particle.model.state)), m)
    end
end


struct Particle{T}
  z::T
end
Particle(z, m) = particleof(m)(z)

function particleof(::AbstractSSMProblem) 

particleof(::LinearSSM) = Vector{Float32}

@yebai
Copy link
Member

yebai commented Aug 3, 2023

@devmotion can you take a look at this PR? I'm very involved in the discussions so might be slightly biased.

Copy link
Member

@yebai yebai left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks, @FredericWantiez -- it looks pretty good to me. I left a lot of minor comments below. Most comments are typos-fixes or clarifications.

In addition, I think we should consider changing.

# Define the structure of the latent space
particleof(::LinearSSM) = Float64
dimension(::LinearSSM) = 2

into something like for clarity and convenience

SSMProblems.particle_eltype(::AbstractSSMModel) = typeof(0.) 
SSMProblems.particle_dimension(::LinearSSM) = 2

So the users only need optionally overload particle_eltype if needed.

EDIT: after introducing the AbstractStateSpaceModel type, I think the AbstractParticle type is redundant and should be removed for maximal simplicity. It is only used for dispatching, which is unnecessary due to the new model argument.

EDIT2: particle_eltype and particle_dimention can be removed too. Such information can be obtained by calling

particle = transition!!(rng, 1, m::LinearSSM, ) # step = 1
particle_eltype = eltype(particle)
particle_dimension = length(particle)

EDIT3: introduce a new transition!! function without the step argument, so it is clear it is for initialization.

src/SSMProblems.jl Outdated Show resolved Hide resolved
src/SSMProblems.jl Outdated Show resolved Hide resolved
src/SSMProblems.jl Outdated Show resolved Hide resolved
src/SSMProblems.jl Outdated Show resolved Hide resolved
src/SSMProblems.jl Outdated Show resolved Hide resolved
docs/src/index.md Show resolved Hide resolved
docs/src/index.md Outdated Show resolved Hide resolved
... # Sample from the transition density
end

function emission_logdensity(step, model::Model, particle::AbstractParticle)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
function emission_logdensity(step, model::Model, particle::AbstractParticle)
function emission_logdensity(step, model::Model, particle::AbstractParticle, observations)

docs/src/index.md Outdated Show resolved Hide resolved
struct Particle{T<:AbstractStateSpaceModel} <: AbstractParticle{T} end

while !all(map(particle -> isdone(t, model, particles), particles)):
ancestors = resample(rng, logweigths)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We need to reset all log weights to zero after the resampling step.

Suggested change
ancestors = resample(rng, logweigths)
ancestors = resample(rng, logweigths)
logweights = zero(logweights)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is resample, also a part of the interface? I think more explanations would make this example much clearer.

src/SSMProblems.jl Outdated Show resolved Hide resolved
docs/src/index.md Outdated Show resolved Hide resolved
src/SSMProblems.jl Outdated Show resolved Hide resolved
docs/src/index.md Outdated Show resolved Hide resolved
Comment on lines +21 to +23
- __Initialisation__: ``f_0(x)``
- __Transition__: ``f(x)``
- __Emission__: ``g(x)``
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it can be a bit confusing here to use x which otherwise only denotes (internal) states. Maybe just define it in/via the model equations below and remove this list?

docs/src/index.md Outdated Show resolved Hide resolved
docs/src/index.md Outdated Show resolved Hide resolved
docs/src/index.md Outdated Show resolved Hide resolved
examples/smc.jl Outdated
end
return trace
return trace[1:(end - 1)] # Discard the root node
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just don't add this node in the first place? The indexing operation here actually creates a copy.

examples/smc.jl Outdated
function linearize(particle::Particle{T,V}) where {T,V}
trace = V[]
current = particle
while !isnothing(current)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In older Julia versions,

Suggested change
while !isnothing(current)
while current !== nothing

was more performant.

examples/smc.jl Outdated
end

ParticleContainer = AbstractVector{<:Particle}
function isdone(t, model::AbstractStateSpaceModel, particles::ParticleContainer)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have my doubts about this function, as mentioned above.

examples/smc.jl Outdated
end

ParticleContainer = AbstractVector{<:Particle}
function isdone(t, model::AbstractStateSpaceModel, particles::ParticleContainer)
return all(map(particle -> isdone(t, model, particle), particles))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Without intermediate allocations1:

Suggested change
return all(map(particle -> isdone(t, model, particle), particles))
return all(particle -> isdone(t, model, particle), particles)

examples/smc.jl Outdated
return logpdf(g(t, particle.state), observations[t])
end

isdone(t, ::LinearSSM, ::AbstractParticle) = t > T
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is an example of a function that closes over a non-const global. Often this leads to very bad performance.

abstract type AbstractParticleCache end

"""
transition!!(rng, step, particle[, cache])
transition!!(rng, model[, timestep, state, cache])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

timestep can be made optional for homogenous SSMs.

# Particle Filter implementation
struct Particle{T} # Here we just need a tree
parent::Union{Particle,Nothing}
# Particle Filter
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we want to keep this code block in the module?


Simulate the particle for the next time step from the forward dynamics.
"""
function transition!! end

"""
transition_logdensity(step, particle, x[, cache])
transition_logdensity(model, timestep, prev_state, next_state[, cache])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
transition_logdensity(model, timestep, prev_state, next_state[, cache])
transition_logdensity(model, timestep, prev_state, current_state[, cache])

src/SSMProblems.jl Outdated Show resolved Hide resolved
Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
@@ -0,0 +1,32 @@
# Concrete Particles as LinkedLists
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@FredericWantiez can we wrap this file into a sub module, e.g. Utils.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Created one Utils, not sure that's the best name

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It follows my lazy naming method. Please feel free to change it if you find a better name.

@yebai yebai merged commit 4153da3 into main Sep 13, 2023
3 checks passed
@yebai yebai deleted the fred/doc branch September 13, 2023 22:30
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants